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ABSTRACT 
The  general  solution  of  the  behavior  of  a  viscously  damped  vibrating 
system  having  several  degrees  of  freedom,  as  developed  in  the  literature 
via  matrix  methods,  is  treated  in  detail  here  and  made  the  basis  of  a 
digital  computer  program  which  is  capable  of  determining  the  natural 
frequencies,  the  mode  shapes,  and  the  displacements  of  each  mass  as 
functions  of  time.  This  program,  which  is  written  in  FORTRAN  60  language 
for  the  Control  Data  Corporation  1604  computer,  affords  several  output 
options.   It  does  not  treat  cases  of  supercritical  damping  or  cases  in 
which  two  or  more  natural  frequencies  are  the  same. 
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Nomenclature 

A  independent  variable  in  reduced  system 

A'  real  part  of  a  complex  matrix 

B  independent  variable  in  reduced  system 

B'  imaginary  part  of  a  complex  matrix 

C  damping  matrix 

D  driving  force  amplitude 

DI  imaginary  component  of  driving  force  amplitude 

DR  real  component  of  driving  force  amplitude 

F(t)  forcing  function 

G  column  of  constant  coefficients 

I  identity  matrix 

K  stiffness  matrix 

M  mass  or  inertia  matrix 

N  degrees  of  freedom 

0  null  matrix 

p  eigenvalue 

R  independent  variable  in  reduced  system 

r,s  subscripts  designating  different  modes 

SS  steady  state  solution 

TS  transient  solution 

U  eigenvector 

V  real  component  of  eigenvector 

W  imaginary  component  of  eigenvector 

X  generalized  displacement,  a  column  vector 

X  generalized  velocity 

iv 


X  generalized  acceleration 

Y  dependent  variable  in  reduced  system 

Z  dependent  variable  in  reduced  system 

o<  decrement  factor   ( 5oon) 

j8  damped  natural  frequency 

§  dependent  variable  in  reduced  system 

Ui  excitation  frequency 

coh  natural  frequency  of  rth  mode 


CHAPTER   I 
INTRODUCTION 
1.1  General  Remarks.   The  analysis  of  a  subcritically  damped  multiple 
degree  of  freedom  vibrating  system  necessitates  obtaining  the  solution  of 
a  complex  eigenvalue  problem  to  determine  the  natural  frequencies  and  mode 
shapes  of  the  system.  Although  the  analysis  presented  in  the  literature 
for  systems  with  two  degrees  of  freedom  may  be  extended  to  systems  with 
more  than  two  degrees  of  freedom,  manual  calculations  are  too  laborious 
to  be  practical.   Therefore  the  natural  frequencies  are  usually  found  by 
ignoring  the  presence  of  the  damping  and  solving  the  resulting  real  eigen- 
value problem.   This  simplifies  the  problem  considerably  and  provides  a 
good  approximation  provided  the  damping  is  light.  Another  technique 
employed  is  to  solve  the  real  eigenvalue  problem  and  then  account  for 
the  damping  using  perturbation  techniques.   However,  even  in  the  absence 
of  damping,  systems  involving  more  than  two  degrees  of  freedom  usually 
require  iteration  or  trial  and  error  techniques  (such  as  the  Stodola  or 
Holzer  methods)  to  obtain  the  mode  shapes. 

The  advent  of  the  electronic  digital  computer  has  eliminated  the 
necessity  of  ignoring  the  damping  component  and  increased  the  size  of  the 
system  for  which  solutions  can  be  obtained.   Although  not  entirely  devoid 
of  error  the  digital  computer  is  highly  reliable  and  its  speed  of  opera- 
tion  has  made  it  an  invaluable  tool  in  engineering  analysis  and  design. 


1.2  Method  of  Analysis.   The  multiple  degree  of  freedom  damped  vibrating 
system  is  described  by  a  set  of  N  linear  second  order  differential  equa- 
tions, where  N  denotes  the  degrees  of  freedom  involved.   Utilizing  matrix 
analysis  and  generalized  coordinates  the  system  is  described  by 

which  is  the  same  form  as  the  single  degree  of  freedom  system.   However, 

•  *• 
M,  C,  and  K  now  represent  square  matrices  and  X,  X,  X,  and  F(t)  are 

column  vectors.   The  analysis  presented  in  what  follows  is  performed  using 
matrix  technique  because  of  the  compactness  of  the  notation  and  the  order- 
ed computational  procedure  which  is  ideally  suited  for  digital  computer 
programming. 

The  mathematical  analysis  is  demonstrated  and  equations  are  derived 
which  describe  the  time  behavior  of  free  and  forced  vibrating  systems. 
Finally  a  digital  computer  program  is  presented  which  performs  the  opera- 
tions indicated  in  the  mathematical  analysis  to  determine  the  natural 
frequencies,  mode  shapes,  and  time  behavior  of  the  damped  vibrating  system. 


CHAPTER  II 


MATHEMATICAL  ANALYSIS 


2.1  Eigenvalues  and  Eigenvectors.  The  viscously  damped  vibrating  system 
with  one  degree  of  freedom  is  described  by  a  linear,  second  order  differ- 
ential equation. 

MX  +  CX  t  KX^  FCO 
M,  C,  and  K  represent  the  mass,  damping  and  stiffness  of  the  system,  and 

F(t)  the  forcing  function.   A  system  with  many  degrees  of  freedom  is  de- 
scribed by  a  set  of  second  order  differential  equations  similar  to  the 
case  of  the  single  degree  of  freedom  system.   Using  matrix  notation  the 
system  is  described  for  the  case  of  free  vibration  by 

(1)      MX  t  C*  +  KX  =  0 

where  M,  C,  K,  are  now  square  matrices  of  order  N,  representing,  as  in 
the  single  degree  case,  the  mass,  damping  and  stiffness  matrices  of  the 
system.   It  is  always  possible  to  write  the  equations  so  that  these  mat- 

•  •  • 

rices  are  symmetric.  N  is  the  number  of  degrees  of  freedom.  X,  X,  and  X 
are  column  vectors  of  order  N  representing  displacement,  velocity,  and 

acceleration  in  generalized  coordinates. 

-1 
Premultiply  equation  (1)  by  M  (Capital  letters  will  hereafter  re- 
present matrices  and  lower  case  letters  scalar  constants)  to  obtain; 

(la)     X  +  M'CX  +  rfV*   =  O 

* 
W.  J.  Duncan  and  A.  R.  Collar  \_l]   have  shown  a  simplified  method  of  find- 
ing the  eigenvalues  and  eigenvectors  by  writing  the  second  order  differen- 
tial equation  in  reduced  form  as  a  first  order  differential  equation.   The 


Numbers  in  brackets  refer  to  references  listed  in  the  bibliography 
page  23 


reduced  form  is  obtained  by  introducing  velocity  as  a  dependent  variable. 

Let 

<2a)      Y- 

Here  Y  is  a  partioned  column  matrix  of  order  2N  with  the  N  displacement 

components  in  the  upper  half  and  N  velocity  components  in  the  lower  half, 

Let 

(2b)      A  ~- 


O 


-i. 


Here  0  is  an  nth-order  null  matrix  and  I  is  an  nth-order  identity  matrix, 
The  reduced  equation  then  becomes 


(3) 


Y  =  AY 


A  solution  of  a  first  order  differential  equation  such  as  equation 
(3)  may  be  taken  as 


(4) 


Y«  u 


Substitution  of  the  assumed  solution  into  equation  (3)  gives 


(5) 


9t  ?t 


Dividing  out  the  exponential  term  and  rearranging,  we  have 


(5a)    JA  -  Pl]  u  -  O 


To  have  a  nontrivial  solution,  it  then  follows  that 

det  |A-pl]*0 
This  is  the  characteristic  matrix  of  the  reduced  system.   The  eigen- 
values or  roots  of  the  characteristic  equation  are  p  ;  r  =  1,  2,  . . . . 2N 


and  the  eigenvectors  are  U  ;  r  =  1,  2,  . ...2N 


where 


u*  = 


*V 


\uz»,r\ 


The  first  subscript  indicates  the  position  in  the  column  and  the  second 
subscript  corresponds  to  the  rtn  eigenvalue. 

As  long  as  the  damping  is  less  than  critical  in  each  mode,  the  re- 
duced system  will  yield  N  complex  conjugate  pairs  of  eigenvalues  and 
eigenvectors.   The  eigenvalues  will  be  of  the  same  form  as  the  sub- 
critically  damped  single  degree  of  freedom  system. 
pr  =  -  Sr  uir    f  j  to,.  /'-.£* 

ft.  *  -Srtor   -  jwr/i-£* 
The  bar  indicates  the  complex  conjugate  and  the  subscript  indicates  the 

mode,  and  where 

Sy.         -  damping  ratio 

to?         -  natural   frequency 

r.  *z 

coryi-*r      -  damped  natural  frequency 
******     -  decrement  factor 
The  corresponding  eigenvectors  are 

2.2  Homogeneous  Equation. 

2.2.1     Orthogonality  Relations.   In  the  case  of  the  undamped  vibrating 
system  it  is  always  possible  to  choose  a  set  of  coordinates  in  which  the 
mass  and  stiffness  matrices  are  uncoupled.   However  in  the  damped  vibrat- 
ing system,  unless  the  damping  matrix  is  proportional  to  either  the  mass 
or  stiffness  matrix,  a  set  of  coordinates  which  will  uncouple  the  equations 
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of  motion  cannot  be  found  without  a  knowledge  of  the  eigenvalues  and  eigen- 
vectors of  the  system. 

K.  A.  Foss  has  developed  a  set  of  orthogonality  relations  for  the 
eigenvectors  of  a  damped  linear  dynamic  system  from  which  coordinates 
which  uncouple  the  equations  of  motion  may  be  found.  |_2J  The  technique  of 
reducing  the  second  order  differential  equation  to  one  of  the  first  order 
is  again  employed,  with  slightly  different  notation. 
Let 


(6a) 


R  = 


O  M 
M  C 


(6  b)   B 


=[- 


-M  O 


K 
Using  the  above  notation,  equation  (1)  becomes 


(6c)    Z  = 


(7) 


TRZ    +■  BH    =  O 


Furthermore   let 


IV* 


Then 


(8) 


where 


x  =  ure 


X     =    PrUr<2 


Z   =  $re 


V* 


fc" 


Ppu, 
uu 


Substituting  equation  (8)  into  equation  (7)  and  dividing  out  the  exponen- 
tial factor  we  get, 


(9) 


PrR$P  +  B$r  =  o 


Since  the  eigenvalues  and  eigenvectors  occur  in  complex  conjugate  pairs 


the  complex  conjugates  of  U  and  p  also  satisfy  the  homogeneous  equation, 


Denoting  the  complex  conjugate  by  the  subscript  S,  equation  (9)  may  be 
written 


(10) 


PsR$>s  +  B$s  s° 


Or  since  M,  C,  and  K  are  symmetric,  R  and  B  must  be  symmetric,  and  equa- 
tion (10)  may  be  written 

(10a)   p5$sTR  f  $38  =  O 

Premultiply  equation  (9)  by   $s   ,  postmultiply  equation  (10a)  by  Q^ 
and  subtract  (10a)  from  (9). 


Cp,-ps)$;r^ 


=  o 


Therefore  unless  r  =  s  the  orthogonality  relations  are 

(id    is*|r  =0 

(12)   §jB$r  *  o 

2.2.2.    General  solution  of  the  homogeneous  equation.   In  the  general 
solution  of  equation  (7),  2N  arbitrary  constants  of  integration  must  be 
evaluated.  The  2N  initial  conditions  of  displacement  and  velocity  are 
used  to  evaluate  these  constants.  Assume  a  solution  of  the  form 


(13)     Em.-^iC.e1 


where  C  represents  the  2N  constants  to  be  determined. 


The  initial  conditions  may  be  expressed  as  a  vector  expansion  of  the 
natural  modes 

(14)  Z(o)=J^     $rC, 

Premultiply  both  sides  of  equation   (14)   by      (g    R 

2  hi 

(15)  $sTRZ(o)  =  2]CR^Cr 


r* 


Using  the  orthogonality  relation,  equation  (11),  the  summation  simplifies 
and  equation  (15)  becomes 


(16)   $^RZCo)  =  £r§sC 


Or 


U7>    C.=  ^l 


The  general  solution  is  therefore  found  by  substituting  the  value  of  the 
constants  of  integration  into  equation  (13). 


(18)    E  tf)  = 


_VAER_2rte* 


§;*$r 


•= i 


It  follows  from  an  evaluation  of  equation  (18)  at  time  zero,  that 

r=  i 

Equation  (18)  may  be  written  in  the  following  partitioned  form  using  the 
nth-order  system  notation. 


(19) 


r.   -i  Y"1 


m?\ 


?rUPUTrMX(o)  f  pJuhuJWXC4  +  prU^UJCXCflS 


from  which 


2M 


(20)  X(*)=, 


V~      [ur.ujM*W>  1-  pr.UrujKAKCo)  *  UruJcX^       ?rV 


rst 


8 


ft 


It  should  be  noted  again  at  this  point  that  for  the  sub-critically  damped 

system  being  considered,  p  and  U  are  complex  quantities  and  occur  in 

r      r 

complex  conjugate  pairs.   The  matrix  operations  indicated  in  equation  (20) 
have  been  carried  out  by  breaking  up  the  complex  quantities  into  real  and 
imaginary  parts  and  then  recombining.  The  operations  are  relatively 
straightforward  but  lengthy,  therefore  only  the  final  form  will  be  pre- 
sented. 

N 

(2i)     x&yjfojF  [^Mxcov+(-^G>rM  +  /srHrM+s»rc)><^cos/erM  e 


hVl^._j  r^MxCo^C-oCt-^H-^GirM  vHrC)  X(a)j  SINA-* 


y--i 


Where: 

ar=  -  Wr(vrTMV,-WrTHWP)-4/3P^TKW,  +VrTCV^  -  W,-TCW,. 

br*      ^/3rCv>Vr-W,TMWr)-4^rVrTMWr    t  £  Vr  C  Wf 
Ar«    VrVj    -WrWrT 

G,.  =  OrAr  +brBr 


The  factor  of  two  in  the  above  equation  results  from  the  reduction  of  the 
summation  from  2N  to  N.   When  the  indicated  operations  of  equation  (20) 
are  carried  out  from  r  =  1  to  r  ■  2N  the  imaginary  components  sum  to  zero; 
therefore  all  quantities  in  equation  (21)  are  real.   Equation  (21)  has 
been  given  by  S.  H.  Crandall  and  R.  B.  McCalley  in  reference  3. 
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2.3  Non-Homogeneous  Equation 

2.3.1     Steady  State  Solution  of  the  Non-Homogeneous  Equation.   Once 

again  using  the  reduced  form,  the  non -homogeneous  equation  becomes 

(22)  RE     *"©2.  «  FW 

Where 

Ftt>« 


o 
In  a  linear  nth-order  differential  equation  the  solution  of  the  non- 


homogeneous  equation  may  be  found,  once  the  solution  of  the  homogeneous 

equation  is  known,  by  expanding  Z      in  a  modal  series.  \2»41 

Thus ,  iti 

(23)    2^  =  ^$^^ 

here   Y,-(fc)    is  a  variable  parameter. 

Upon  substitution  of  equation  (23)  in  equation  (22),  we  obtain 

(24)     y-RbX  +  ^§ry,-  =  F^ 

-T 
Premultiplication  of  equation  (24)  by   $r     and  making  use  of  the 

orthogonality  relations,  yields 

(25)     $Tp*$r*.  +  $>$ryr «  $>« 

Or 

(25a)        Rny    -  ?rRoy  -    Fn 

Where 

Fn  *      C  F  C\) 

10 


Using  the  Laplace  Transform  technique  the  solution  of  equation  (25a) 

is  easily  determined. 

Rearranging 

(a)   Y-p.Y  .-^- 

The  Laplace  Transform  is 

(b>    syci)  -  ?rvu)  =  £v£S 

Therefore 

(O   Y<s)  -  ^  (5_M 

The  inverse  transform  of  (c)  is 

(26)  Yto.J-.y*e*C    F,(r)dt 

o 
Since    pr  =  -<Xr  +j  £r 

(26a)    Yft)»  -^-J  £*"  *"*  pnC^)  (cosprfr-tj+j  6ii4|8rft-t))at 

O 

Substitute  (26a)  into  (23) 

(27)      Z($  =  )2&j£*r     '    £(*)  (cosfrQt-*)  ¥  jSiM/Srfr-t^dt 


t-=i 


As  in  the  case  of  the  homogeneous  solution  the  reduced  system  quantities 
are  replaced  by  their  nth-order  equivalence.   The  final  form  of  the 


steady  state  solution  then  becomes 


(28)  X^^>/^)C  ^"^  C**+C**Wt 


r 


/ qj+bf 


-tfr(W) 


+  \  Zik.    fCt)  e~  SH4j8r(^-tWC 


The  notation  in  equation  (28)  is  the  same  as  that  used  in  the  homogeneous 
solution. 
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2.3.2     Steady  State  Solution  with  Sinusoidal  Excitation.   A  special, 
but  very  common  problem  is  the  case  in  which  the  forcing  function,  f  £V)   , 
is  sinusoidal.   Expressing  the  forcing  function  as 


(29)   fctt  =  tflie^) 


where    GL   -  indicates  "the  real  part  of" 

D  -  column  of  driving  force  amplitudes,  possibly  complex  to  admit 
phasing 

OJ  -  excitation  frequency 

The  non-homogeneous  equation  therefore  becomes 

(30)         MX     +  CX    f  KX    »    C3t(l5CJt0  ) 
Assume  a  solution  of  the  form 


(31) 


X(VJ  -  £  (GejUJ  ) 


where  G       is  a  column  of  undetermined  constant  coefficients.   Substitute 
the  assumed  solution  in  equation  (30)  and  divide  through  by  the  exponen- 
tial factor  to  obtain 

(32)  -o/-MG,    +jcoCG,    +  KG     -"B 

Or 

-l 

(33)  G    =  [k-co^M  +j  uC]l 

Therefore  the  steady  state  solution  for  the  case  of  sinusoidal  excitation 
is 


(34)    XW=  (K.[K-^MtjcuC]T)CJ 
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2.4  Transient  Plus  Steady  State  Solution. 

2.4.1.   Initial  Conditions.   The  total  solution  is  given  by 

(35)  Xto  <=  T5  tSS 

where  TS  is  the  so  called  transient  solution  and  SS  the  so  called  steady 

state  solution.   However  in  the  total  solution  the  initial  conditions 

of  displacement  and  velocity  in  the  transient  solution  must  account  for 

the  initial  values  in  the  steady  state  solution.  This  is  accomplished 

as  follows.   In  the  case  of  sinusoidal  excitation  the  steady  state 

solution  was  found  to  be 

.i 


(34)     X(r)=   &(    K-o>**M  +  J6oC]DeJUi 


To  allow  for  phase  differences  in  the  exciting  forces  on  the  masses,  let 
D  be  given  by 

(36)  B  =  DR  +jDI 

and  let  the  inverse  of  IK-  co^M-^JcjCJ  be  denoted  by 

(37)  A     +  j  "&' 

Then  since    £       -   CoScot    ¥  j  SisJcot  equation   (34) 

becomes 

(38)  Xto  ^  ^J[a+jb][^^X5i](co5co^+JS^^^J 
Carrying  out   the  indicated  multiplications,    we  obtain 

(39)  X  (t)  =  (  (a'dr-BBI)  cos  tor  -  (b'DR  +  At>l)  sin  cot  J 
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and  the  velocity  is  given  by 

(40)  X  (V)  -j-6o(ADR-B'Dl)  siNcor  -  co  (b'DR  4- a'di)  Cos  co  ^  J 

Therefore  the  initial  values  of  displacement  and  velocity  of  the  steady 
state  solution  are  obtained  by  evaluating  equations  (38)  and  (39)  at  ta  ° 

(41)  Xss(o)  =  A'DR  -o'dI 

(42)  *sSCo)  =  -co  (b'dK.  +  a'di) 

The  initial  conditions  used  in  determining  the  time  behavior  of  the 
transient  portion  for  the  total  solution  therefore  become 

(43)  Xc  Co)  ^  *(g)  -  (a'D*.  "  Sbl) 

(44)  vcCo)  -  XCo)  +■  6o  (Vt>R  4-  a'di) 

The  complete  time-solution  is  thus  obtained  as  indicated  in  equation  (35) 
where  the  transient  solution  is  found  by  using  the  modified  initial  condi- 
tions of  equations  (43)  and  (44). 
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CHAPTER  III 
PROGRAM  DESCRIPTION 
3.1  General  Remarks.   A  digital  computer  program  "PROGRAM  VIBSYS"  is 
presented  which  performs  the  mathematical  operations  of  the  equations 
developed  in  the  previous  chapter.   The  program  is  coded  in  Fortran  60 
programming  language,  ^5, 6J  specifically  for  the  Control  Data  Corpora- 
tion 1604  computer.  Although  the  Fortran  60  language  is  applicable  to 
most  large  digital  computers  there  are  minor  variations  which  are  pecul- 
iar to  the  specific  system  in  use.  The  Fortran  60  language  does  not  per- 
mit automatic  operations  with  complex  numbers;  therefore  all  operations 
involving  complex  numbers  are  accomplished  by  operating  on  the  real  and 
imaginary  parts  separately  and  then  recombining. 

The  eigenvalues  and  eigenvectors  of  the  reduced  system  are  found  by 
a  matrix  iteration  scheme  which  utilizes  the  direct  and  inverse  power 
methods  and  matrix  deflation.   A  mathematical  subroutine  MATSUB  is  used 
to  carry  out  these  operations.  [7,81  As  presented,  a  maximum  of  twenty 
iterations  will  be  performed  using  the  direct  power  method  and  then  a 
maximum  of  twenty  using  the  inverse  power  method.   If  convergence  is  not 
reached  with  the  inverse  power  method  a  print  out  to  this  effect  will  be 
executed. 

Subroutine  "INVERT"  is  used  for  matrix  inversion.    "INVERT"  uses 
the  Gaussian  elimination  and  pivotal  techniques.   Inversion  of  the 
M K.-u/'M  +  .'coCJ   matrix  which  contains  complex  elements  is  achieved  in 


Subroutine  INVERT  is  a  library  subroutine  of  the  Naval  Postgraduate 
School  and  is  designated  locally  as  Fl  NPS  INVERT. 
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the   following  manner.        Let    \k  +  jB]    be   the  matrix  to  be   inverted  and 
\C  +  jDj    the  inverse   to  be  determined.      Then  by  definition 

[A+j6][c+jT5]=  I 
The  above  matrix  equation  leads   to   two  simultaneous   equations  with 

two  unknowns,    C  and  D. 

(a)  AC-BD  =    I 

(b)  -QC  +AD  =    O 
Solving   first   for  C 

c=    [a  +  Wb] 

and   then  for  D 

d  =  -  a'bc 

Thus  the  complex  inversion  is  reduced  to  multiplications,  additions, 
and  inversions,  of  matrices  of  real  numbers. 

Flexibility  and  utility  are  the  principal  aims  of  the  program. 
Usage  requires  a  knowledge  of  the  mass,  stiffness,  and  damping  matrices. 
Although  various  output  options  are  available  which  require  additional 
input  data,  the  three  above  mentioned  matrices  are  all  that  are  necessary 
to  determine  the  eigenvalues  and  eigenvectors  of  the  reduced  system  and 
the  natural  frequencies  and  mode  shapes  of  the  original  system. 

3.2  Program  Options.  In  addition  to  finding  the  natural  frequencies  and 
mode  shapes  of  the  system,  five  options  are  available  which  describe  the 
time  behavior  of  the  system  under  conditions  of  free  and  forced  vibration. 
(a)  Option  1.   The  time  solution  of  the  free  vibration  problem 
is  obtained  in  general  form  and  no  additional  input  data  is  required  for 
execution  of  this  option.   In  section  2.2.2  the  general  solution  of  the 
homogeneous  equation  was  developed  and  is  repeated  here  for  convenience. 
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r-=.l 
The  output  of  option  1  consists  of  the  four  coefficient  matrices  of  the 
XtoCoSfSr^j  Xto)Cosfr):    )   X  C«>  SiNi^t  and  X<^  SvH^rt 

terms.   Therefore,  the  output  of  option  1  will  consist  of  4N  square  matrices. 

(b)  Option  2.   The  execution  of  option  2  requires  as  additional 
input  data  the  values  of  the  initial  displacement,  and  initial  velocity 
vectors.   The  product  of  the  coefficient  matrices  of  option  1  and  the 
initial  displacement  vector  or  initial  velocity  vector,  as  appropriate,  is 
performed  to  obtain  a  coefficient  column  for  the  cosine  and  sine  terms. 

(c)  Option  4.*  The  general  steady  state  solution  of  the  forced 
vibration  problem  is  provided  by  option  4.   The  general  solution  with  a 
forcing  function  f(t)  was  developed  in  section  2.3.1  and  found  to  be 

(28)   xft)»)  J^r-  JfCE)e""   "rcoS/3ra-t)<=lT 


I-.1 


kJ  * 


z 


r-s  i 
In  this  case  the  coefficient  matrices  of  the  convolution  integral  are 

evaluated.   There  will  be  2N  square  matrices,  N  corresponding  to  the 

coefficient  of  the  integral  involving  the  cosine  term  and  N  coefficient 

matrices  of  the  integral  involving  the  sine  term. 

(d)   Option  5.   The  steady  state  solution  of  the  special  case 

of  forced  vibration  with  sinusoidal  excitation  is  determined  in  option  5. 


* 
Option  3  is  described  in  subsection  (e),  following. 
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In  section  2.4.1  the  expression  for  the  time  behavior  was  found  to  be 

(39)      x(\)  -    (ADR-BDl)  cos co V    -  (T3DRtADl)  si^tot 

Additional  input  data  required  for  the  execution  of  this  option  include 
the  driving  force  amplitude  and  excitation  frequency.   The  output  consists 
of  the  coefficient  column  vectors  of  the  cosine  and  sine  terms. 

(e)  Option  3.   A  plot  of  displacement  versus  time  is  made  for 
each  mass,  for  one  of  the  three  following  cases. 

1.  Free  Vibration 

2.  Steady  state  vibration  with  sinusoidal  excitation 

3.  Transient  plus  steady  state  with  sinusoidal  excita- 
tion 

Case  1  is  obtained  when  the  initial  conditions  of  displacement  and  veloc- 
ity are  given  as  input  data  and  the  driving  force  amplitudes  and  excita- 
tion frequency  are  zero* 

Case  2  is  obtained  when  the  driving  force  amplitudes  and  excitation 
frequency  are  given  and  the  initial  conditions  of  displacement  and 
velocity  are  zero. 

Case  3  is  obtained  when  the  initial  conditions,  driving  force  amplitudes 
and  excitation  frequency  are  all  given. 

Any  combination  of  the  options  may  be  obtained  with  one  set  of  input 
data.   Input  data  format  and  output  control  is  described  in  detail  in 
Appendix  B. 

3.3  Accuracy  of  Method.   The  accuracy  of  the  solution  is  contingent  on 
the  accuracy  with  which  the  eigenvalues  and  eigenvectors  are  determined 
and  the  effect  of  computer  roundoff  error.   Internal  checks  have  been 
provided  so  that  the  integrity  of  the  solution  may  readily  be  evaluated. 
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3.3.1  Internal  Checks.   Prior  to  the  determination  of  the  eigenvalues 
and  eigenvectors  the  trace  and  determinate  of  the  matrix,  A,  is  calcu- 
lated.  A  check  is  therefore  readily  available  since  the  trace  of  the 
matrix  A  is  equal  to  the  sum  of  the  eigenvalues  of  the  characteristic 
matrix  of  A  and  the  determinant  is  equal  to  the  product  of  the  eigen- 
values.  The  trace  and  determinant  of  A  and  the  sum  and  product  of  the 
eigenvalues  are  included  as  standard  output  in  Program  VIBSYS. 

Additional  checks  are  inherent  in  the  solution.   In  Option  1  the 
summation  of  the  coefficient  matrices  of  the  X(o)  Cos/3Ki  terms 

must  sum  to  the  identity  matrix,  since  the  evaluation  of  x(t)  at  time 
zero  must  equal  the  initial  conditions  of  displacement.   This  is  easily 
seen  by  referring  to  equation  21  and  considering  the  initial  values  of 
velocity  to  be  zero.   Similarly,  in  Option  2  the  sum  of  the  coefficient 
columns  of  the  cosine  terms  must  equal  the  initial  conditions  of  dis- 
placement. 

In  the  steady  state  solution  with  sinusoidal  excitation, the  inversion 
of  the  complex  matrix  presents  one  of  the  greatest  possible  sources  of 
error.  Therefore  when  Option  5  is  executed,  the  output  includes  the  real 
and  imaginary  parts  of  the  inverse  and  the  product  of  the  original  matrix 
and  its  inverse.   This  product  should,  of  course,  be  the  identity  matrix. 

3.3.2  External  Checks.   In  order  to  determine  the  reliability  of 
Program  VIBSYS,  sample  problems  of  two  degrees  of  freedom  for  which  solu- 
tions were  available  in  the  literature  were  programmed.   Correlation  of 
results  was  excellent.   No  suitable  sample  problems  involving  more  than 
two  degrees  of  freedom  were  found  in  the  literature.   Therefore,  in  order 
to  extend  the  range  of  testing,  solutions  of  undamped  systems  were  compar- 
ed with  program  solutions  of  similar  systems  with  negligible  damping. 

19 


The  results  were  as  anticipated;  as  the  damping  was  decreased  the  natural 
frequencies  approached  those  of  the  undamped  system.   The  results  of  a 
sample  problem  are  demonstrated  in  Appendix  E. 

3.4   Program  Limitations.   Although  the  digital  computer  extends  the  size 
of  system  for  which  solutions  are  obtainable,  the  ultimate  size  is  dictat- 
ed by  the  storage  capacity  of  the  computer.   The  CDC  1604  computer  has 
approximately  32,000  storage  locations;  however,  it  is  easily  seen  that 
this  may  be  rapidly  exhausted.   Program  VIBSYS  is  limited  to  a  system 
with  10  degrees  of  freedom. 

The  basic  principle  in  the  development  of  the  general  equations  is 
the  orthogonality  relations  of  the  eigenvectors  and  the  parameters  of 
the  system.  Therefore  each  eigenvalue  must  have  associated  with  it  a 
unique  eigenvector.  This  requirement  is  not  fulfilled  in  the  case  of 
a  system  having  two  equal  eigenvalues,  and  therefore  two  natural  fre- 
quencies of  equal  magnitude.  In  such  a  case  VIBSYS  determines  the 
frequencies  and  the  program  terminates  at  that  point. 

If  the  damping  is  critical  in  any  part  of  the  system  an  attempt 
will  be  made  to  obtain  the  eigenvalues  and  eigenvectors.   However,  if 
they  are  found,  they  will  not  occur  in  the  complex  conjugate  pairs  and 
the  program  will  be  terminated. 
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CHAPTER  IV 
CONCLUSIONS  AND  RECOMMENDATIONS 

Program  VIBSYS  provides  an  accurate  method  for  determining  the  natural 
frequencies,  mode  shapes,  and  time  behavior  of  a  subcritically  damped 
dynamic  system.   The  compiling  time  for  the  program  is  4  minutes  and  30 
seconds,  while  the  running  time  for  a  system  with  four  degrees  of  freedom 
is  approximately  one  minute.   (Running  time  will  depend  on  the  number  of 
iterations  required  to  obtain  the  eigenvalues  and  eigenvectors.)  There- 
fore it  is  more  efficient,  with  respect  to  computer  time,  to  make  multiple 
runs. 

Comparison  of  natural  frequencies  of  undamped  and  lightly  damped 
systems  has  shown  that  the  undamped  approximation  is  very  good  for  damping 
ratios  of  less  than  0.01.   For  systems  tested,  up  to  four  degrees  of 
freedom,  the  natural  frequencies  of  the  damped  and  undamped  systems  did 
not  differ  in  the  first  three  significant  figures. 

The  program  may  be  extended  in  a  variety  of  possible  ways.   The 
general  nature  of  Program  VIBSYS  requires  liberal  use  of  computer  stor- 
age space.   The  size  of  the  system  may  be  extended  by  segmenting  the 
program  into  separate  programs  for  handling  specific  problems  such  as, 
free  vibration,  forced  vibration  with  a  general  forcing  function,  or 
forced  vibration  with  sinusoidal  excitation.   Furthermore,  with  minor 
modifications  the  present  program  may  be  terminated  upon  finding  the 
natural  frequencies  and  mode  shapes  and  used  as  a  separate  program  to 
obtain  input  data  for  specialized  programs. 

Program  VIBSYS  can  be  augmented  and  made  more  useful  by  eliminating 
the  restrictions  which  cause  the  program  to  be  terminated  if  either  (a) 
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the  system  has  two  eigenvalues  of  equal  magnitude,  or  (b)  one  of  the  modes 
has  aperiodic  free  motion,  that  is  a  purely  real  eigenvalue. 
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APPENDIX  A 
PROGRAM  STRUCTURE 
A. 1  General  Remarks.   Program  VIBSYS  is  composed  of  a  main  body  which  is 
divided  into  five  sections,  and  two  subroutines,  INVERT  and  MATSUB.   Func- 
tion subroutines  used  include  the  sine,  cosine,  exponential,  square  root, 
and  absolute  value  functions.   These  function  subroutines  are  called  from 
the  Fortran  library  tape.   Subroutine  draw  is  also  available  on  the  For- 
tran library  type  and  is  programmed  specifically  for  the  CDC  1604  computer 
for  use  with  the  CALCOMP  plotter. 

Subroutine  MATSUB,  which  is  available  as  a  CO-OP  Library  mathematical 
subroutine,  required  minor  modifications  in  order  to  retain  the  eigen- 
values and  eigenvectors  in  storage  for  use  in  the  main  body  of  the  program. 
The  flow  chart  provided  with  MATSUB  in  the  CO-OP  Library  was  found  to  be 
inadequate  for  understanding.   In  order  to  perform  the  necessary  modifica- 
tions a  revised  flow  chart  was  made  up  and  included  in  Appendix  C. 

A.  2  Main  Body.   The  main  body  consists  of  five  sections  with  functions 
as  listed  below 

(a)  Input  -  In  addition  to  allocation  of  storage  spaces  for 
dimensioned  arrays  and  setting  up  a  control  for  reading  input  data,  the 
input  section  contains  a  control  sequence  for  MATSUB  and  curve  labeling. 

(b)  Natural  frequencies  and  mode  shapes  -  The  reduced  character- 
istic matrix  is  constructed  and  the  eigenvalues  and  eigenvectors  are  deter- 
mined.  The  natural  frequencies  and  mode  shapes  are  then  found. 

(c)  Forced  Vibration  -  The  coefficient  column  vectors  of  the 
sine  and  cosine  terms  of  option  5  are  calculated. 

(d)  Free  Vibration  -  The  outputs  of  options  1,  2,  and  3  are 
formulated. 
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(e)  Forced  Vibration  with  General  Excitation  -  The  coefficient 
matrices  of  option  4  are  calculated. 

A. 3   Subroutines. 

(a)  MATSUB  -  Evaluates  the  eigenvalues  and  eigenvectors  of  the 
reduced  system. 

(b)  INVERT  -  Inverts  a  real  matrix. 
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AB 

ABC 

ALIS 

ALRS 

AM 

AMD 

AMM 

AMS 

ARN 

ARY 

BRN 

CD 

CDC 

CPD 

DI 

DISPI 

DISPR 

DM 

DR 

EP1 

EP2 


PROGRAM  NOTATION 
coefficient  column  vector  of  cos/8r  t  term,  output  of 
option  2 

coefficient  column  vector  of  cos/6r  t  term  when  X(0),  and 
X(0)  are  modified  to  obtain  total  solution 
control  parameter  for  MATSUB 
control  parameter  for  MATSUB 
inverse  of  mass  matrix 

product  of  inverse  of  mass  matrix  and  damping  matrix 
mass  matrix  lb-sec  /in 

product  of  inverse  of  mass  matrix  and  stiffness  matrix 
-2  o^(VrTMVr  -  WrTMWr)  -4  (Br  J  VrTMW.  +  VrTCVr  -  WrTCWr 
imaginary  part  of  identity  matrix 
2  /9r  (VrTMVr  -  WrTMWr)  -  4  c»<rVrTMWr  +  2VrTCVr 
coefficient  column  vector  of  sin  /St   term,  output  of 
option  2 

coefficient  column  vector  of  si.n/8   t  term  when  X(0)  and 
X(0)  are  modified  to  obtain  total  solution 
coefficient  column  of  sincot,  output  of  option  5 
imaginary  part  of  driving  force  amplitude 
imaginary  part  of  modal  matrix   ,* 
real  part  of  modal  matrix 
damping  matrix  lb°sec/in 
real  part  of  driving  force  amplitude 
iteration  parameter  for  MATSUB 
iteration  parameter  for  MATSUB 
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GBI 

GBR 

GRM 

HRM 

IC 

IDET 

IEG 

IFC 

ITITLE 

IVEC 

MIT 

MITS 

MP1 

MP2 

MP3 

MP4 

MP5 

N 

RIN 

RMI 

RPI 

S 

SPEC 

STEP 

TIPI 


iteration  parameter  for  MAT SUB 
iteration  parameter  for  MATSUB 

coefficient  matrix  of  X(0)  cos/9   t,  output  of  option  1 
coefficient  matrix  of  X(0)  sin/Qrt,  output  of  option  1 
control  parameter  for  VIBSYS 
control  parameter  for  MATSUB 
control  parameter  for  MATSUB 
control  parameter  for  VIBSYS 
title  for  graphical  output 
control  parameter  for  MATSUB 

controls  number  of  iterations  in  MATSUB,  power  method 
controls  number  of  iterations  in  MATSUB  inverse  power 
method 

option  1  control 
option  2  control 
option  3  control 
option  4  control 
option  5  control 
order  of  system 
real  part  of  identity  matrix 

coefficient  column  vector  of  cos  co  t  term,  output  of 
option  5 

real  part  of  the  inverse  of   [k-  coTl  +  j  to  6\ 
stiffness  matrix  lb/in 
spectral  matrix 

step  size  for  graphical  output 
imaginary  part  of  inverse  of   K- co  m  +  j  o  C 
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UI 

UR 

VALI 

VALR 

VCO 

VCV 

VCW 

VECI 

VECR 

VO 

VMV 

VMW 

WCW 

WDM 

WMK 

WMW 

WO 

X 

XCO 


matrix  of  reduced  system,  imaginary  part 

matrix  of  reduced  system,  real  part 

imaginary  part  of  eigenvalue  of  reduced  system  (/S^) 

real  part  of  eigenvalue  of  reduced  system    (  °^) 

initial  velocity  vector  modified  to  account  for  initial 

velocity  of  steady  state  solution 

vrTcvr 
vrTcwr 

imaginary  part  of  eigenvector  of  reduced  system  (V  ) 

real  part  of  eigenvector  of  reduced  system  (W  ) 

initial  velocity  vector 

VrTMVr 

VrTMWr 

wrTcwr 

imaginary  part  of  \K   -  OJ T4  +  j  60  Cj 
real  part  of   [K  -  6J  m  +  j6iCj 
WrTMWr 


excitation  frequency 

time  coordinate  for  graphs 

initial  condition  of  displacement  modified  to  account  for 

initial  displacement  of  steady  state  solution 
XO        -    initial  displacement  vector 

XOC       -    coefficient  matrix  of  X(0)  cos/Srt  term,  output  of  option  1 
XOS       -    coefficient  matrix  of  X(0)  sin^j-t  term,  output  of  option  2 
Y        -    ordinate  of  graphical  output 

The  above  notation  lists  the  principal  array  and  parameter  names  used  in 
the  main  body  of  the  program.   Array  names  not  listed  are  used  for  inter- 
mediate operations.  Where  appropriate  the  names  are  associated  with  the 
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symbols  used  in  the  mathematical  analysis, 
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APPENDIX  B 


INSTRUCTIONS  FOR  USE  OF  PROGRAM 


B.l  Purpose.   The  purpose  of  Program  VIBSYS  is  to  determine  the  natural 
frequencies,  mode  shapes,  and  time  behavior  of  a  subcritically  damped, 
linear  dynamic  system  with  N  degrees  of  freedom.   (N  Z.    10)   The  mass, 
damping  and  stiffness  matrices  of  the  system  must  be  known. 
B.2   Input  Data.   A  blank  card  must  follow  the  last  end  card  of  the  pro- 
gram deck.   The  data  cards  then  follow  in  the  order  and  format  listed  be- 
low.  The  first  data  card  uses  the  input  format  915  and  is  the  control 
card  for  the  program.   The  nine  fields  are  designated  as  follows: 


1. 
2. 


N 

MP1  = 


the  order  of  the  system 


3.   MP2 


4.   MP3 


5.   MP4  = 


6.   MP5 


7. 
8. 
9. 


IC 


IFC  = 


NPTS 


Option  1  executed 
Omit 

Option  2  executed 
Omit 

Option  3  executed 
Omit 

Option  4  executed 
Omit 

Option  5  executed 
Omit 

Initial  conditions  given 
Initial  conditions  are  zero 

Amplitude  and  frequency  of  excitation  given 
Omit 


Number  of  points  to  be  calculated  for  graphical 
output 


All  nine  fields  must  be  right  justified. 
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The  control  card  is  followed  by  the  mass,  damping,  and  stiffness 
matrices  respectively.   Each  matrix  is  read  in  rowwise  using  input  format 
8F10.3*  with  each  row  starting  on  a  new  card  for  systems  with  N  ^  8.   For 
systems  with  N  >  8  the  elements  are  read  in  rowwise  in  a  sequential  manner, 
that  is  the  first  element  of  the  second  row  should  appear  in  the  field  ad- 
jacent to  the  last  element  of  the  first  row.   In  this  case  each  matrix  must 
start  on  a  new  card. 

The  initial  displacement  vector,  initial  velocity  vector,  real  part 
of  the  driving  force  amplitude  vector,  and  imaginary  part  of  the  driving 
force  amplitude  vector,  follow,  in  that  order,  using  format  8F10.3.* 
Each  vector  must  start  on  a  new  card. 

The  final  data  card  uses  a  2F10.3*  format.   The  first  field  contains 
the  excitation  frequency  and  the  second  the  step  size  (i.e.,  time  increment) 
for  the  graphical  output  option. 

Blank  cards  must  be  used  for  initial  conditions  of  displacement  and 
velocity,  real  and  imaginary  parts  of  forcing  function  amplitudes  and 
final  data  card  when  these  values  are  zero.  Multiple  runs  may  be  process- 
ed by  placing  additional  data  decks  behind  the  first.   The  program  is 
terminated  by  placing  a  blank  card  behind  the  last  data  card. 


* 
Should  this  format  present  undesirable  restrictions  on  the  size  of 
the  elements  of  the  matrix  the  so  called  "E"  field  may  be  used  by  changing 
card  number  32.   The  "E"  field  was  avoided  since  it  is  more  susceptible  to 
error  by  the  user. 
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B.3  Deck  Assembly.   The  first  card  of  the  program  deck  is  a  job 
card.   The  statement  "Use  scratch  tape"  must  be  included  in  addition  to 
the  required  job  card  information. 

JOB  CARD 

PROGRAM  VIBSYS 


END 

"subroutine  INVERT 


END 

SUBROUTINE  MATSUB 


END 

END 

BLANK  CARD 

DATA  CARDS 


BLANK  CARD 
B.4  Cautions  to  User.  The  curves  drawn  in  the  graphical  output  option 
are  straight  line  approximations  between  computed  values.   The  step  size 
must  therefore  be  chosen  appropriately  in  order  to  obtain  a  smooth  curve. 
Since  the  maximum  number  of  points  permitted  by  subroutine  DRAW  is  900, it 
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may  not  be  possible  to  see  the  transient  phase  completely  die  out. 

If  difficulty  is  encountered  in  obtaining  the  eigenvalues  and  eigen- 
vectors, the  control  parameters  of  MATSUB  may  be  altered  to  enhance  con- 
vergence.  The  value  of  MIT  and  MITS,  card  numbers  20  and  21,  control  the 
maximum  number  of  iterations  to  be  performed  in  the  power  method  and  in- 
verse power  method  respectively.   ALRS  and  ALIS,  card  numbers  78  and  79, 
represent  a  complex  parameter  that  may  be  called  the  "origin".   MATSUB 
will  usually  converge  on  the  eigenvalue  most  distant  from  the  "origin". 
ALRS  =  1.0  and  ALIS  =0.1  in  Program  VIBSYS.   The  convergence  criteria 
for  the  power  method  is  set  by  the  value  of  EP1,  card  number  85  while 
the  inverse  power  method  is  determined  by  EP2  card  number  86.   EP1  ■  10~\ 
and  EP2  =  10~14  in  the  present  version.   Furthermore,  if  IEG,  card  number 
76,  is  set  equal  to  one,  the  eigenvalue  iterants  will  be  printed  out  and 
more  suitable  values  of  ALRS  and  ALIS  may  possibly  be  found  by  inspection. 
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APPENDIX  C 


FLOW  CHARTS 


The  symbols  used  in  the  flow  charts  are  defined  as  follows: 


o 


arithmetic,  read  or  print  operations 


connectors 


-  decision 


call  subroutine  indicated 


The  notation  used  in  the  flow  chart  for  MATSUB  corresponds  to  that  used 
in  the  CO  OP  Library  version,  reference  7. 
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FLOW  CHART     MAM  BODY    PROGRAM  VIBSYS 


START 


READ:  PROGRAM  PARAMETERS.  N,  MP].,  MP2, 
MP3.  MP4,  MP5.  IC.  IFC,  &PTS 


READ:  MASS,  DAMPING,  AND 
.  STIFFNESS  MATRICES 


PRINT:  MASS,  DAMPING,  AND 
STIFFNESS  MATRICES 


READ:  INITIAL  CONDITIONS 
XO#  VO 


PRINT:  XO,  YO 


I 


READ:  DRIVING  FORCE  AMPL. 
DR,  DI 


PRINT:  DR,  DI 


READ:  EXCIT.  FREQ.  AND 
STEP  SIZE.  WO,  STEP 


PRINT:  WO,  STEP 


I 


READ:  ITERATION  PARAMETERS 


J_ 


READTCURVE-  LABELS 


END  OF  INPUT  DATA 
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14 


AMS-  -AM*S 


45- 


AMD=-AM*DM 


FORM  UR  MATRIX 


I 


UI=0.0 


I  CALL  MATSUB 

\uivui 


PRINT:  TRACE,  DET,  EIGENVALUES, 
EIGENVECTORS,  SUM,  AND  PROD 


SELECT  POSITIVE  VALUES  OF  VALI 


i 


SPEC=VALI**2-VALR**Z 


31 


VJNAT  =  SQRTFCll  SPEC II    ) 


SELECT  VECR,  VECI  CORRESPONDING  TO 
POSITIVE  VALI 


I 


DISPR  -  VECR 
DISPI =  VECI 


PRINT:  NATURAL  FREQUENCIES 
AND  MODE  SHAPES 
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YES 


951 


PRINT  WNAT  =  0  SYSTEM  IS 
PRORART.Y  OVF.RDAMPP.D 


3£ 


GO  TO 


lnoo. 


KO 


START  FORCED  VIBRATION  WITH 
SINUSOIDAL  EXCITATION 


NO 


DIFs  WNAT(I)-  WNAT(J) 


tfigf  S-WO**2*AMM 


501 


WDM  =  WO*  DM 


YES 


941 


PRINT  SISTEH  HAS 
FKEQ  OF  EQUAL  MAGN  . 


GO  TO  1000 


WMXI  =  WMK 


k^CALL  INVERT^ 
\  503 


WIC  =  WMKI*WDM 


|          504 

cwic  =  wdm*wic 

■ 

RPI  =  WMK  +  CWIC 

" 

TIPI=  -WIC*RPI 

...         ii 

CHECK  INVERSION  OF  K  -  w2K  t  IwC 

' 

807 

R1N  =  A1U  -  BID 
ARY  =  ETC  -  AID 

510 

SMI  =  RTD  -  TIDI 
CPD  =  RTD  +  TTD 

v 

START  FREE  VIBRATION  PROBLEM 
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I  -  1,M 


VMV  =  DI3PR*AMM*DISPR 
WMW  =  DISPI*AMM*DISPI 
VMW"  =  DISPK*AMM*DISPI 
WCW  =  DI3PI*DM*DISPI 
VCV  =  DISPE*DM*DISPB 
VCW  =  DISPR*DM*DISPI 


JZ1 


CALCULATE  ARNfH?W 


80 


81 


GR  =  b/ARN**2  +  BRN**2 
KR  =  -HM/ARH**2  +  BRH**2 


I 


6RH  =  2*GR*AMM 
HRM  =  2*HR*AMM 
GRC  a 2*GR*DM 
HRC  =  2*HR*DM 


I 


82 


XOC  =  -  AGM  +  £HM  +  GRC 
X03-  -BGM  -AHM  +  HRC 


NO 

_<^F 

' 

■ 

201 

AB  =  GRI«l*VO  +  XOC*XO 
CD  =  HRM*VO  *■  X03*X0 

=  1? 


YES 

■     101 

! 

KC  =  1 

i 

' 

PRINT     XCC,XOS 

. 

1     ' 

GO  TO  100 

i 

MP1  =  0 

■ 

. 

GO  TO  60 
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<C^P3  ~  lj^> 

NO 

1 

1    YES 

GO  TO  400 

KC=  3 

'( 

709 

xco  =  xo  -  rke 

VCO  =  VO  -  WO*CPD 

" 

710 

ABC  =  GRM*VCO  +  XCC*XCO 

cdc  -  hrm*vco  +  xos*xco 

J  711 

AsGR  i'ABC 


CcHIIIfllEL 


I=1,N 


NO 


YES 


insL 


CALCULATE   GRAPH  POINTS  FOR 
T]?AMsti.-mt  SOLUTION 


CALCULATE   GRAPH  POINTS  FOR 
STEADY  STATE  SOLUTION 
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CALL  DRAW 


410 


GR  = 
HR  = 

2* 
2* 

(G 

*B)/ARN**2  +  EKN**2 

-  HM)  /ARN**2  +  ifiN**2 

V 

PRINT 

GR,HK 

NO 

<^Q?5  *l?^ 

>_X1^ 

1 

f 

PRINT  RPI,TIPI,RIN, 
ARY.RM-CPD 

'' 

~ 

GO  XO  1000 

NO       ^ 

■ 

t 

READ  NEXT  SET  OF 
DATA  CARDS 

C-7 


SUBROUTINE  MATSUB 


SET:  M,  IEG,  XVEC,  <*„  , 
/9  ,   MIT,  MITS,<?(  ,  £2  ' 

■ 

' 

SUM  =  0 
PROD  =  0 

' 

COMPUTE  TRACE  OF  C 

■ 

• 

C-A 
C—B 

■ 

■   . 

ASSIGN  520  TO  IA 
ASSIGN  811TO  ID 
INTER  0 

TRIANGULARIZE  A 


COMPUTE  DET  A 


_1 


PRINT  TRACER,TRACEI, 
DET  A 


ASSIGN  912  TO  ID 
ASSIGN  530  TO  IA 
ASSIGN  40  TO  IB 
ASSIGN  523  TO  IC 


B-A 

"T 

ISLs  0 


SET  <* 

I 
IT  -  1 

i 

X  =   1.0 


(A-«*I)-»A 


C-8 


1  =  1 

5 


ITERATION  POWER  METHOD 


YES 


SET  *'UvO 
ICI-  1QQQ 


PRINT  N,ALR,ALI,ICT,IJ 


X;* NORMALIZED  Y; 


VALR(N)  =  ALR 
YALI(N)^  ALI 


ICT  =  i 

MIT2  =  HITS  +  i 


<*i-  «*«,+■  m; 


d 


C-9 


© 


INVERSE  PQV/Eft  METHOD 
I 310 


(A-^I)-^A 


I 


U=IJ+  1 


2$. 


V 


TRIANGULARIZE  A 


YES 


i^-m, 


I 


SET  FLAG 
tht=2000 

\ 


SET  UP  FOR  EACK 
SUBSTITUTION 


JL 


SMaLL=  1000 


_Z52 


Y:    (A-^i)y-  =  x 


K-\ 


\  =  RnSMETZED  Y/ 


535 


530 


NO 


1  40 


IZZ=  1 


-42- 


^2_ 


PRINT  2,1,*, 


32 


B->  A 


RAYLEIGH  QUOTIENT 


YES 


AX^  ,  xa 


_50_ 


NO 


C-10 


YES 
_fiJL 


ISL  =  0 


X 


PRINT  N,*:  ,   ICT,   i 


SAVE    <*; 
SUM  =  SUM  +  *> 
PROD  5  PROD*  g<j 


X 


N  *  N  - 1 


MATRIX  DEFLATION 


m  -  N 


B=  DEFLATED  A 


(NO 

402 

<Ki  = 

o<-  + 

e 

1 

IT- 

I-*-! 

I 

PRINT    ©< 

. 

YES 


(C-^I)->A 


INTER  =  0 

ASSIGN  013  TO  TA 


YES  .fR 


PRINT  SUM, PROD 


T 


990 


_£XIT_ 


67 


*i  =  °H,1 


SUM  =  SUM  +  <x; 
PROD-  PROD* *£ 
SAVR  *t 


"     523 


ISL=  0 


$ 


PRINT 


©<:, 


VALR(l)  -ALR 
,VAU(1)=ALI 


\     320 


N=  0 

NILL.1 


5 
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TRIANOULASIZE  A 


INTER ^  1 
ASSIGN  530  TO  IA 


NO 


535 


911 


YES 


520 


COMPUTE  DET  A 


912 


PRINT  PET  A 


924 


COMPUTE  AND  PRINT  LARGEST 
AND  SMALLEST  DIAGONAL 

ELEHKNTS 


— >;«« 


J 


YES 


X=0.0 


ASSIGN  753  TO  IB 
ASSIGN  704  TO  IC 

♦   31  - 


YR  1.0 

YI  ,Q,Q 


914 


ISL=  -1 


IVEC  = 1? 


BACKWARD  SUBSTITUTION 


704 


PRINT  EIGENVECTOR 


I 


VECR(I,NN)=XR(I) 
VECI(I.NN)=XI(I) 


ASSIGN  40  TO  IB 


NO 


ASSIGN  525  TO  IC 


B-^A 


¥ 
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APPENDIX  D 


PROGRAM  LISTING 
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APPENDIX  E 


SAMPLE  PROBLEM 


The  sample  problem  presented  is  that  of  a  four  degree  of  freedom 
system  shown  in  the  figure  below.   In  order  to  show  a  comparison  of  data 
for  damped  and  undamped  systems,  the  example  chosen  is  one  in  which  the 
natural  frequencies  of  the  undamped  system  may  be  obtained  with  relative 
ease.   The  program  is  not  restricted  to  problems  with  the  symmetry  dis- 
played in  the  example. 


Let       ML  =  M2  =  M3  -  M4  =  10  lb  sec/in 

Kx  ■  K2  -  K3  =  K4  =  1000  lb/in 

Cj_  =  C2  =  C_  =  C4  =  30  lb  sec/in 
For  the  undamped  case  the  frequency  determinate  then  becomes 


-   KU>  *■ 

-K 

0 

0 

-K 

2K  -  MdO 

-K 

0 

0 

-K 

2K  -  MO>2 

-K 

0 

0 

-K 

2K 

-    M0l> 

E-l 


Expanding  the  determinate 
6 


K\  6 
The  roots  of  the  frequency  equation  are 


K£  4. 


CO,  =  6. 18034, £0^  =  11. 75571, ^3  =  16. 18034, 604=  19.02113 

The  problem  was  programmed  with  a  forcing  function  f (t)  =  lOOsin 
(60t) ,  and  an  initial  displacement  on  M4of  0.5  inches.   A  step  size  of 
0.005  seconds  was  used  for  the  graphical  option.   All  five  options  were 
called  for  and  the  resulting  output  is  shown  on  the  following  pages. 

In  addition  to  the  output  presented  for  C  =  30  lbsec/in  runs  were  made 
with  C  =  0.01,  0.10,  and  1.0  lb  sec/in.   The  frequencies  and  damping 
ratios  for  all  runs  are  tabulated  below. 


C  lb  sec/in 

$, 

CO,  ,  rad/sec 

3 
^2. 

Cd2  ,rad/sec 

0.01 

— 

6.18034 

0.0001 

11.75571 

0.10 

0.0003 

6.18034 

0.0006 

11.75570 

1.00 

0.0031 

6.18028 

0.0067 

11.75530 

30.00 

0.0098 

6.12699 

0.1821 

11.38430 

C 

lb  sec/in 

33 

^  ^  rad/sec 

0.0001 

£0   rad/sec 

0.01 

0.0001 

16.18034 

19.02113 

0.10 

0.0008 

16.18033 
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